#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Thomas Breuer.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file contains methods for general linear mappings of finite
##  dimensional free left modules.
##
##  There are two default representations of such general mappings,
##  one by generators and images, the other by two bases and a matrix.
##
##  Note that a matrix is not the appropriate object to represent a general
##  linear mapping if it is not total or not single-valued;
##  moreover, if one does not prescribe images of a basis but of an
##  arbitrary generating system, one does not want to compute a basis at the
##  time the general mapping is formed;
##  finally, the matrix is not appropriate to compute preimages, whereas
##  the general mapping by images behaves symmetrically in this respect.
##
##  (The matrix is best for linear mappings used as arithmetic elements,
##  for mapping elements of the source to the range and back; storing
##  images and preimages avoids the matrix multiplication.)
##
##  1. methods for linear general mappings given by images
##  2. methods for linear mappings given by matrices
##  3. methods for vector spaces of linear mappings
##  4. methods for algebras of linear mappings
##  5. methods for full hom spaces
##

#T TODO:
#T
#T specific representation for nat. hom.
#T     (choice of coefficients instead of silly matrix)
#T AsLeftModuleGeneralMappingByImages, to allow id + id, c * id, -id,
#T     Zero( id ), id + zero, - zero, c * zero, ...
#T \= methods for m.b.m. and g.m.b.i. (if bases coincide, compare data)
#T parent dependencies for nat. hom.

#T put bases into mappings;
#T check that they are really bases of source/range!


#############################################################################
##
##  1. methods for linear general mappings given by images
##

#############################################################################
##
#R  IsLinearGeneralMappingByImagesDefaultRep
##
##  is a default representation of $F$-linear general mappings between two
##  free left modules $V$ and $W$ where $F$ is equal to the left acting
##  domain of $V$ and of $W$.
##
#T  (It would be possible to allow situations where $F$ is only contained
#T  in the left acting domain of $W$;
#T  this would lead to asymmetry w.r.t. taking the inverse general mapping.)
##
##  Defining components are
##
##  `generators' \: \\
##      list of vectors in $V$,
##
##  `genimages' \: \\
##      list of vectors in $W$.
##
##  The general mapping is defined as the linear closure of the relation
##  that joins the $i$-th entry in `generators' and the $i$-th entry in
##  `genimages'.
##
##  If one wants to compute images, one needs the components
##  `basispreimage' \: \\
##      a basis of the $F$-module generated by `generators',
##
##  `imagesbasispreimage' \: \\
##      images of the basis vectors of `basispreimage',
##
##  `corelations' \: \\
##      linearly independent generators for the corelation space,
##      i.e., of the space of all row vectors <r> such that
##      `LinearCombination( <r>, generators )' is zero in $V$.
##      (The corresponding linear combinations of `genimages'
##      generate the cokernel.)
##
##  If these components are not yet bound, they are computed by
##  `MakeImagesInfoLinearGeneralMappingByImages' when they are needed.
##  If `generators' is a *basis* of a free left module then these
##  components can be entered without extra work.
##
##  If one wants to compute preimages, one needs the components
##  `basisimage' \: \\
##      a basis of the $F$-module generated by `genimages',
##
##  `preimagesbasisimage' \: \\
##      preimages of the basis vectors of `basisimage',
##
##  `relations' \: \\
##      linearly independent generators for the relation space,
##      i.e., of the space of all row vectors <r> such that
##      `LinearCombination( <r>, genimages )' is zero in $W$.
##      (The corresponding linear combinations of `generators'
##      generate the kernel.)
##
##  If these components are not yet bound, they are computed by
##  `MakePreImagesInfoLinearGeneralMappingByImages' when they are needed.
##  If `genimages' is a *basis* of a free left module then these
##  components can be entered without extra work.
##
##  Computed images and preimages of free left modules under linear mappings
##  are always free left modules.
##  If one needs more structure (e.g., that of an algebra) for an image or
##  preimage then the linear mapping must have a special representation.
##
##  Note that the inverse general mapping of a linear mapping defined by
##  images is best handled if it uses the default method,
##  since such an inverse general mapping delegates the tasks of computing
##  (pre)images to the original general mapping.
##  So the (pre)images info is computed only once.
#T  but what about sums of such mappings?
#T  better try to share info also in this case?
#T  (share a list that is filled with the info later?)
##
DeclareRepresentation( "IsLinearGeneralMappingByImagesDefaultRep",
    IsAttributeStoringRep,
    [ "basisimage", "preimagesbasisimage", "corelations",
      "basispreimage", "imagesbasispreimage", "relations",
      "generators", "genimages" ] );

InstallTrueMethod( IsAdditiveElementWithInverse,
    IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep );
InstallTrueMethod( IsLeftModuleGeneralMapping,
    IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep );


#############################################################################
##
#M  LeftModuleGeneralMappingByImages( <S>, <R>, <gens>, <imgs> )
##
InstallMethod( LeftModuleGeneralMappingByImages,
    "for two free left modules and two homogeneous lists",
    [ IsFreeLeftModule, IsFreeLeftModule,
      IsHomogeneousList, IsHomogeneousList ],
    function( S, R, gens, imgs )

    local map;        # general mapping from <S> to <R>, result

    # Check the arguments.
    if   Length( gens ) <> Length( imgs )  then
      Error( "<gens> and <imgs> must have the same length" );
    elif not IsSubset( S, gens ) then
      Error( "<gens> must lie in <S>" );
    elif not IsSubset( R, imgs ) then
      Error( "<imgs> must lie in <R>" );
    elif LeftActingDomain( S ) <> LeftActingDomain( R ) then
      Error( "<S> and <R> must have same left acting domain" );
    fi;

    # Make the general mapping.
    map:= Objectify( TypeOfDefaultGeneralMapping( S, R,
                             IsSPGeneralMapping
                         and IsLeftModuleGeneralMapping
                         and IsLinearGeneralMappingByImagesDefaultRep ),
                     rec() );

    SetMappingGeneratorsImages(map,[gens,imgs]);
    # Handle the case that `gens' is a basis.
    if IsBasis( gens ) then
      map!.basispreimage       := gens;
      map!.imagesbasispreimage := imgs;
      map!.corelations         := Immutable( [] );
    fi;

    # Handle the case that `imgs' is a basis.
    if IsBasis( imgs ) then
      map!.basisimage          := imgs;
      map!.preimagesbasisimage := gens;
      map!.relations           := Immutable( [] );
    fi;

    # return the general mapping
    return map;
    end );


#############################################################################
##
#M  LeftModuleHomomorphismByImagesNC( <S>, <R>, <gens>, <imgs> )
##
InstallMethod( LeftModuleHomomorphismByImagesNC,
    "for two left modules and two lists",
    [ IsFreeLeftModule, IsFreeLeftModule, IsList, IsList ],
    function( S, R, gens, imgs )
    local map;        # homomorphism from <source> to <range>, result
    map:= LeftModuleGeneralMappingByImages( S, R, gens, imgs );
    SetIsSingleValued( map, true );
    SetIsTotal( map, true );
    return map;
    end );


#############################################################################
##
#F  LeftModuleHomomorphismByImages( <S>, <R>, <gens>, <imgs> )
##
InstallGlobalFunction( LeftModuleHomomorphismByImages,
    function( S, R, gens, imgs )
    local hom;
    hom:= LeftModuleGeneralMappingByImages( S, R, gens, imgs );
    if IsMapping( hom ) then
      return LeftModuleHomomorphismByImagesNC( S, R, gens, imgs );
    else
      return fail;
    fi;
end );


#############################################################################
##
#M  AsLeftModuleGeneralMappingByImages( <linmap> )  . for a lin. gen. mapping
##
InstallMethod( AsLeftModuleGeneralMappingByImages,
    "for a linear g.m.b.i.",
    [     IsLeftModuleGeneralMapping
      and IsLinearGeneralMappingByImagesDefaultRep ],
    IdFunc );


#############################################################################
##
#M  ImagesSource( <map> ) . . . . . . . . . . . . . . . . for linear g.m.b.i.
##
InstallMethod( ImagesSource,
    "for a linear g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    if IsBound( map!.basisimage ) then
      return UnderlyingLeftModule( map!.basisimage );
    else
      return SubmoduleNC( Range( map ), MappingGeneratorsImages(map)[2] );
#T is it used that the second argument may be a basis object?
    fi;
    end );


#############################################################################
##
#M  PreImagesRange( <map> ) . . . . . . . . . . . . . . . for linear g.m.b.i.
##
InstallMethod( PreImagesRange,
    "for a linear g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    if IsBound( map!.basispreimage ) then
      return UnderlyingLeftModule( map!.basispreimage );
    else
      return SubmoduleNC( Source( map ), MappingGeneratorsImages(map)[1] );
#T is it used that the second argument may be a basis object?
    fi;
    end );


#############################################################################
##
#F  MakeImagesInfoLinearGeneralMappingByImages( <map> )
##
##  Provide the information for computing images, that is, set up
##  the components `basispreimage', `imagesbasispreimage', `corelations'.
##
BindGlobal( "MakeImagesInfoLinearGeneralMappingByImages", function( map )
    local preimage,
          ech,
	  mapi,
          B;

    preimage:= PreImagesRange( map );
    mapi:= MappingGeneratorsImages( map );

    if   Dimension( preimage ) = 0 then

      # Set the entries explicitly.
      map!.basispreimage       := Basis( preimage );
      map!.corelations         := IdentityMat( Length( mapi[2] ),
                                      LeftActingDomain( preimage ) );
      map!.imagesbasispreimage := Immutable( [] );

    elif IsGaussianRowSpace( Source( map ) ) then
#T operation MakeImagesInfo( map, source )
#T to leave this to the method selection ?
#T or flag `IsFromGaussianSpace' ?

      # The images of the basis vectors are obtained on
      # forming the linear combinations of images of generators
      # given by `ech.coeffs'.

      ech:= SemiEchelonMatTransformation( mapi[1] );
      map!.basispreimage       := SemiEchelonBasisNC(
                                      preimage, ech.vectors );
      map!.corelations         := Immutable( ech.relations );
      map!.imagesbasispreimage := Immutable( ech.coeffs * mapi[2] );
#T problem if mapi[2] is a basis and if this does not store that it is a small list!

    else

      # Delegate the work to the associated row space.
      B:= Basis( preimage );
      ech:= SemiEchelonMatTransformation( List( mapi[1],
                     x -> Coefficients( B, x ) ) );
      map!.basispreimage       := BasisNC( preimage,
                                      List( ech.vectors,
                                        x -> LinearCombination( B, x ) ) );
      map!.corelations         := Immutable( ech.relations );
      map!.imagesbasispreimage := Immutable( List( ech.coeffs,
                                        x -> LinearCombination( x,
                                                 mapi[2] ) ) );

    fi;
end );


#############################################################################
##
#F  MakePreImagesInfoLinearGeneralMappingByImages( <map> )
##
##  Provide the information for computing preimages, that is, set up
##  the components `basisimage', `preimagesbasisimage', `relations'.
##
BindGlobal( "MakePreImagesInfoLinearGeneralMappingByImages", function( map )
    local image,
          ech,
	  mapi,
          B;

    mapi:= MappingGeneratorsImages( map );
    image:= ImagesSource( map );

    if   Dimension( image ) = 0 then

      # Set the entries explicitly.
      map!.basisimage          := Basis( image );
      map!.relations           := IdentityMat( Length( mapi[1] ),
                                      LeftActingDomain( image ) );
      map!.preimagesbasisimage := Immutable( [] );

    elif IsGaussianRowSpace( Range( map ) ) then

      # The preimages of the basis vectors are obtained on
      # forming the linear combinations of preimages of genimages
      # given by `ech.coeffs'.
      ech:= SemiEchelonMatTransformation( mapi[2] );
      map!.basisimage          := SemiEchelonBasisNC( image, ech.vectors );
      map!.relations           := Immutable( ech.relations );
      map!.preimagesbasisimage := Immutable( ech.coeffs * mapi[1]);
#T problem if mapi[1] is a basis and if this does not store that it is a small list!

    else

      # Delegate the work to the associated row space.
      B:= Basis( image );
      ech:= SemiEchelonMatTransformation( List( mapi[2],
                     x -> Coefficients( B, x ) ) );
      map!.basisimage          := BasisNC( image,
                                      List( ech.vectors,
                                        x -> LinearCombination( B, x ) ) );
      map!.relations           := Immutable( ech.relations );
      map!.preimagesbasisimage := Immutable( List( ech.coeffs,
                                      row -> LinearCombination(
                                                 row, mapi[1] ) ) );

    fi;
end );


#############################################################################
##
#M  CoKernelOfAdditiveGeneralMapping( <map> ) . . .  for left module g.m.b.i.
##
InstallMethod( CoKernelOfAdditiveGeneralMapping,
    "for left module g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local genimages;

    # Form the linear combinations of the basis vectors for the
    # corelation space with the `genimages' of `map'.

    if not IsBound( map!.corelations ) then
      MakeImagesInfoLinearGeneralMappingByImages( map );
    fi;
    genimages:= MappingGeneratorsImages(map)[2];
    return SubmoduleNC( Range( map ),
               List( map!.corelations,
                     r -> LinearCombination( genimages, r ) ) );
    end );


#############################################################################
##
#M  IsSingleValued( <map> ) . . . . . . . . . . . .  for left module g.m.b.i.
##
InstallMethod( IsSingleValued,
    "for left module g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local genimages;

    if not IsBound( map!.corelations ) then
      MakeImagesInfoLinearGeneralMappingByImages( map );
    fi;
    genimages:= MappingGeneratorsImages(map)[2];
    return ForAll( map!.corelations,
                   r -> IsZero( LinearCombination( genimages, r ) ) );
    end );


#############################################################################
##
#M  KernelOfAdditiveGeneralMapping( <map> ) . . . .  for left module g.m.b.i.
##
InstallMethod( KernelOfAdditiveGeneralMapping,
    "for left module g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local generators;

    # Form the linear combinations of the basis vectors for the
    # relation space with the `generators' of `map'.

    if not IsBound( map!.relations ) then
      MakePreImagesInfoLinearGeneralMappingByImages( map );
    fi;
    generators:= MappingGeneratorsImages(map)[1];
    return SubmoduleNC( Source( map ),
               List( map!.relations,
                     r -> LinearCombination( generators, r ) ) );
    end );


#############################################################################
##
#M  IsInjective( <map> )  . . . . . . . . . . . . .  for left module g.m.b.i.
##
InstallMethod( IsInjective,
    "for left module g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local generators;

    if not IsBound( map!.relations ) then
      MakePreImagesInfoLinearGeneralMappingByImages( map );
    fi;
    generators:= MappingGeneratorsImages(map)[1];
    return ForAll( map!.relations,
                   r -> IsZero( LinearCombination( generators, r ) ) );
    end );


#############################################################################
##
#M  ImagesRepresentative( <map>, <elm> )  . . . . .  for left module g.m.b.i.
##
InstallMethod( ImagesRepresentative,
    "for left module g.m.b.i., and element",
    FamSourceEqFamElm,
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep,
      IsObject ],
    function( map, elm )
    if not IsBound( map!.basispreimage ) then
      MakeImagesInfoLinearGeneralMappingByImages( map );
    fi;
    elm:= Coefficients( map!.basispreimage, elm );
    if elm = fail then
      return fail;
    elif IsEmpty( elm ) then
      return Zero( Range( map ) );
    fi;
    return LinearCombination( map!.imagesbasispreimage, elm );
    end );


#############################################################################
##
#M  PreImagesRepresentative( <map>, <elm> ) . . . .  for left module g.m.b.i.
##
InstallMethod( PreImagesRepresentative,
    "for left module g.m.b.i., and element",
    FamRangeEqFamElm,
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep,
      IsObject ],
    function( map, elm )
    if not IsBound( map!.basisimage ) then
      MakePreImagesInfoLinearGeneralMappingByImages( map );
    fi;
    elm:= Coefficients( map!.basisimage, elm );
    if elm = fail then
      return fail;
    fi;
    return LinearCombination( map!.preimagesbasisimage, elm );
    end );


#############################################################################
##
#M  ViewObj( <map> )  . . . . . . . . . . . . . . .  for left module g.m.b.i.
##
InstallMethod( ViewObj,
    "for a left module g.m.b.i",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local mapi;

    mapi:= MappingGeneratorsImages( map );
    View( mapi[1] );
    Print( " -> " );
    View( mapi[2] );
    end );


#############################################################################
##
#M  PrintObj( <map> ) . . . . . . . . . . . . . . .  for left module g.m.b.i.
##
InstallMethod( PrintObj,
    "for a left module g.m.b.i",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local mapi;

    mapi:= MappingGeneratorsImages( map );
    Print( "LeftModuleGeneralMappingByImages( ",
	   Source( map ), ", ", Range( map ), ", ",
	   mapi[1], ", ", mapi[2], " )" );
    end );

InstallMethod( PrintObj,
    "for a left module hom. b.i",
    [ IsMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local mapi;

    mapi:= MappingGeneratorsImages( map );
    Print( "LeftModuleHomomorphismByImages( ",
	   Source( map ), ", ", Range( map ), ", ",
	   mapi[1], ", ", mapi[2], " )" );
    end );


#############################################################################
##
#M  \*( <c>, <map> )  . . . . . . . . . . . .  for scalar and linear g.m.b.i.
##
InstallMethod( \*,
    "for scalar and linear g.m.b.i.",
    [ IsMultiplicativeElement,
      IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( scalar, map )
    local mult,   # the multiple of `map', result
	  mapi,	  # generators and images
          F;      # left acting domain

    # Check the scalar.
    # (Maybe it is in fact another mapping, and we want to compose.)
    if     not IsInt( scalar )
       and not IsElmsColls( FamilyObj( scalar ),
                   FamilyObj( LeftActingDomain( Range( map ) ) ) ) then
      TryNextMethod();
    fi;

    mapi:=MappingGeneratorsImages(map);
    # Construct the linear general mapping (if possible).
    mult:= LeftModuleGeneralMappingByImages(
               Source( map ), Range( map ), mapi[1],
               List( mapi[2], v -> scalar * v ) );

    # Maintain info on the preimage side of the general mapping.
    if IsBound( map!.basispreimage ) then
      mult!.basispreimage       := map!.basispreimage;
      mult!.imagesbasispreimage := Immutable(
          List( map!.imagesbasispreimage, v -> scalar * v ) );
      mult!.corelations         := map!.corelations;
    fi;

    # Being a mapping is preserved by scalar multiplication.
    if HasIsSingleValued( map ) then
      SetIsSingleValued( mult, IsSingleValued( map ) );
    fi;
    if HasIsTotal( map ) then
      SetIsTotal( mult, IsTotal( map ) );
    fi;

    # If the scalar is invertible in the left acting domain of the source
    # then surjectivity and injectivity are maintained as well as the image.
    F:= LeftActingDomain( Source( map ) );
    if scalar in F and IsUnit( F, scalar ) then

      if HasIsInjective( map ) then
        SetIsInjective( mult, IsInjective( map ) );
      fi;
      if HasIsSurjective( map ) then
        SetIsSurjective( mult, IsSurjective( map ) );
      fi;

      if IsBound( map!.basisimage ) then
        scalar:= Inverse( scalar );
        mult!.basisimage          := map!.basisimage;
        mult!.preimagesbasisimage := Immutable(
            List( map!.preimagesbasisimage, v -> scalar * v ) );
        mult!.relations           := map!.relations;
      fi;

    fi;

    return mult;
    end );


#############################################################################
##
#M  AdditiveInverseOp( <map> )  . . . . . . . . . . . . . for linear g.m.b.i.
##
InstallMethod( AdditiveInverseOp,
    "for linear g.m.b.i.",
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map )
    local ainv,   # the additive inverse of `map', result
          mapi;

    mapi:=MappingGeneratorsImages(map);
    # Construct the linear general mapping (if possible).
    ainv:= LeftModuleGeneralMappingByImages(
               Source( map ), Range( map ), mapi[1],
               List( mapi[2], AdditiveInverse ) );

    # Maintain images and preimages info.
    if IsBound( map!.basispreimage ) then
      ainv!.basispreimage       := map!.basispreimage;
      ainv!.imagesbasispreimage := Immutable(
          List( map!.imagesbasispreimage, AdditiveInverse ) );
      ainv!.corelations         := map!.corelations;
    fi;
    if IsBound( map!.basisimage ) then
      ainv!.basisimage          := map!.basisimage;
      ainv!.preimagesbasisimage := Immutable(
          List( map!.preimagesbasisimage, AdditiveInverse ) );
      ainv!.relations           := map!.relations;
    fi;

    # Being a mapping is preserved by scalar multiplication.
    if HasIsSingleValued( map ) then
      SetIsSingleValued( ainv, IsSingleValued( map ) );
    fi;
    if HasIsTotal( map ) then
      SetIsTotal( ainv, IsTotal( map ) );
    fi;

    # Surjectivity and injectivity are maintained.
    if HasIsInjective( map ) then
      SetIsInjective( ainv, IsInjective( map ) );
    fi;
    if HasIsSurjective( map ) then
      SetIsSurjective( ainv, IsSurjective( map ) );
    fi;

    return ainv;
    end );


#############################################################################
##
#T  \<( <map1>, <map2> )
##
##  method for two linear mappings from Gaussian spaces, use canonical bases?
##


#############################################################################
##
#M  CompositionMapping2( <map2>, map1> )   for left mod. hom. & lin. g.m.b.i.
##
InstallMethod( CompositionMapping2,
    "for left module hom. and linear g.m.b.i.",
    FamSource1EqFamRange2,
    [ IsLeftModuleHomomorphism,
      IsLeftModuleGeneralMapping
      and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map2, map1 )
    local comp,        # composition of <map2> and <map1>, result
	  mapi1,
          gens,
          genimages;

    # Check that the linear mappings can be composed.

    mapi1:=MappingGeneratorsImages(map1);
    # Compute images for the generators of `map1'.
    if     IsLinearGeneralMappingByImagesDefaultRep( map2 )
       and mapi1[2] = MappingGeneratorsImages(map2)[1] then

      gens      := mapi1[1];
      genimages := MappingGeneratorsImages(map2)[2];

    else

      gens:= mapi1[1];
      genimages:= List( mapi1[2],
                        v -> ImagesRepresentative( map2, v ) );

    fi;

    # Construct the linear general mapping.
    comp:= LeftModuleGeneralMappingByImages(
               Source( map1 ), Range( map2 ), gens, genimages );

    # Maintain images info (only if `gens' is not a basis).
    if     IsLinearGeneralMappingByImagesDefaultRep( comp )
       and not IsBound( comp!.basispreimage  )
       and IsBound( map1!.basispreimage ) then
      comp!.basispreimage       := map1!.basispreimage;
      comp!.corelations         := map1!.corelations;
      comp!.imagesbasispreimage := Immutable(
          List( map1!.imagesbasispreimage,
                v -> ImagesRepresentative( map2, v ) ) );
    fi;

    # Return the composition.
    return comp;
    end );


#############################################################################
##
#M  \+( <map1>, map2> ) . . . . . . . . . . . . . . . for two linear g.m.b.i.
##
##  If both general mappings respect zero, additive inverses, scalar
##  multiplication then the sum also does.
##
InstallOtherMethod( \+,
    "for linear g.m.b.i. and general mapping",
    IsIdenticalObj,
    [ IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep,
      IsGeneralMapping ],
    function( map1, map2 )
    local gens,
          genimages,
	  mapi1,
          sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    elif  PreImagesRange( map1 ) <> PreImagesRange( map2 ) then
      Error( "<map1> and <map2> must have same preimage" );
    fi;

    mapi1:=MappingGeneratorsImages(map1);

    if     IsLinearGeneralMappingByImagesDefaultRep( map2 )
       and mapi1[1] = MappingGeneratorsImages(map2)[1] then

      # If the generators in both general mappings are the same,
      # it suffices to add the images.
      gens      := mapi1[1];
      genimages := mapi1[2] + MappingGeneratorsImages(map2)[2];

    else

      # Compute images of the generators of `map1' under `map2'.
      # (Note that both general mappings must be described in terms of
      # `generators' in order to keep the meaning of `corelations'.)
      gens:= mapi1[1];
      genimages:=   mapi1[2]
                  + List( mapi1[1],
                          v -> ImagesRepresentative( map2, v ) );

    fi;

    # Construct the linear general mapping.
    sum:= LeftModuleGeneralMappingByImages(
              Source( map1 ), Range( map1 ), gens, genimages );

    # Maintain images info (only if `gens' is not a basis).
    if     IsLinearGeneralMappingByImagesDefaultRep( sum )
       and IsLinearGeneralMappingByImagesDefaultRep( map2 )
       and not IsBound( sum!.basispreimage  )
       and IsBound( map1!.basispreimage )
       and IsBound( map2!.basispreimage )
       and map1!.basispreimage = map2!.basispreimage then
      sum!.basispreimage       := map1!.basispreimage;
      sum!.corelations         := map1!.corelations;
      sum!.imagesbasispreimage :=
          map1!.imagesbasispreimage + map2!.imagesbasispreimage;
    fi;

    # Return the sum.
    return sum;
end );

InstallOtherMethod( \+,
    "for general mapping and linear g.m.b.i.",
    IsIdenticalObj,
    [ IsGeneralMapping,
      IsGeneralMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map1, map2 )
    local gens,
          genimages,
	  mapi2,
          sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    elif  PreImagesRange( map1 ) <> PreImagesRange( map2 ) then
      Error( "<map1> and <map2> must have same preimage" );
    fi;

    mapi2:=MappingGeneratorsImages(map2);

    if     IsLinearGeneralMappingByImagesDefaultRep( map1 )
       and MappingGeneratorsImages(map1)[1]= mapi2[1] then

      # If the generators in both general mappings are the same,
      # it suffices to add the images.
      gens      := mapi2[1];
      genimages := MappingGeneratorsImages(map1)[2] + mapi2[2];

    else

      # Compute images of the generators of `map1' under `map2'.
      # (Note that both general mappings must be described in terms of
      # `generators' in order to keep the meaning of `corelations'.)
      gens:= mapi2[1];
      genimages:=   List( mapi2[1],
                          v -> ImagesRepresentative( map1, v ) )
                  + mapi2[2];

    fi;

    # Construct the linear general mapping.
    sum:= LeftModuleGeneralMappingByImages(
              Source( map1 ), Range( map1 ), gens, genimages );

    # Maintain images info (only if `gens' is not a basis).
    if     IsLinearGeneralMappingByImagesDefaultRep( sum )
       and IsLinearGeneralMappingByImagesDefaultRep( map1 )
       and not IsBound( sum!.basispreimage  )
       and IsBound( map1!.basispreimage )
       and IsBound( map2!.basispreimage )
       and map1!.basispreimage = map2!.basispreimage then
      sum!.basispreimage       := map1!.basispreimage;
      sum!.corelations         := map1!.corelations;
      sum!.imagesbasispreimage :=
          map1!.imagesbasispreimage + map2!.imagesbasispreimage;
    fi;

    # Return the sum.
    return sum;
end );


#############################################################################
##
#M  \+( <map1>, <map2> )  . . . . . . . . . for two linear mappings by images
##
##  The method for (total and single-valued general) mappings takes
##  advantage from the fact that `generators' and `basispreimage' components
##  need not be distinguished since the `corelations' component is empty.
##
InstallOtherMethod( \+,
    "for linear m.b.i. and mapping",
    IsIdenticalObj,
    [ IsMapping and IsLinearGeneralMappingByImagesDefaultRep,
      IsMapping ],
    function( map1, map2 )
    local gens,
          genimages,
	  mapi1,
          sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    elif  PreImagesRange( map1 ) <> PreImagesRange( map2 ) then
      Error( "<map1> and <map2> must have same preimage" );
    fi;

    if     IsBound( map1!.basispreimage ) then

      # Use the basis in the construction.
      gens:= map1!.basispreimage;

      if     IsLinearGeneralMappingByImagesDefaultRep( map2 )
         and IsBound( map2!.basispreimage )
         and map1!.basispreimage = map2!.basispreimage then

        genimages := map1!.imagesbasispreimage + map2!.imagesbasispreimage;

      else

        genimages:=   map1!.imagesbasispreimage
                    + List( gens,
                            v -> ImagesRepresentative( map2, v ) );

      fi;

    else

      mapi1:=MappingGeneratorsImages(map1);

      if     IsLinearGeneralMappingByImagesDefaultRep( map2 )
	  and mapi1[1] = MappingGeneratorsImages(map2)[1] then

	# If the generators in both general mappings are the same,
	# it suffices to add the images.
	gens      := mapi1[1];
	genimages := mapi1[2] + MappingGeneratorsImages(map2)[2];

      else

	# Compute images of the generators of `map1' under `map2'.
	# (Note that both general mappings must be described in terms of
	# `generators' in order to keep the meaning of `corelations'.)
	gens:= mapi1[1];
	genimages:=   mapi1[2]
		    + List( mapi1[1],
			    v -> ImagesRepresentative( map2, v ) );

      fi;
    fi;

    # Construct the linear mapping.
    sum:= LeftModuleHomomorphismByImagesNC(
              Source( map1 ), Range( map1 ), gens, genimages );

    # Return the sum.
    return sum;
    end );

InstallOtherMethod( \+,
    "for mapping and linear m.b.i.",
    IsIdenticalObj,
    [ IsMapping,
      IsMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    function( map1, map2 )
    local gens,
          genimages,
	  mapi1,mapi2,
          sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    elif  PreImagesRange( map1 ) <> PreImagesRange( map2 ) then
      Error( "<map1> and <map2> must have same preimage" );
    fi;

    if     IsBound( map2!.basispreimage ) then

      # Use the basis in the construction.
      gens:= map2!.basispreimage;

      if     IsLinearGeneralMappingByImagesDefaultRep( map1 )
         and IsBound( map1!.basispreimage )
         and map1!.basispreimage = map2!.basispreimage then

        genimages := map1!.imagesbasispreimage + map2!.imagesbasispreimage;

      else

        genimages:=   List( gens, v -> ImagesRepresentative( map1, v ) )
                    + map2!.imagesbasispreimage;

      fi;

    else

      mapi2:=MappingGeneratorsImages(map2);

      if     IsLinearGeneralMappingByImagesDefaultRep( map1 )
	  and MappingGeneratorsImages(map1)[1] = mapi2[1] then

	# If the generators in both general mappings are the same,
	# it suffices to add the images.
	gens      := mapi2[1];
	genimages := MappingGeneratorsImages(map1)[2] + mapi2[2];

      else

	# Compute images of the generators of `map2' under `map1'.
	# (Note that both general mappings must be described in terms of
	# `generators' in order to keep the meaning of `corelations'.)
	gens:= mapi2[1];
	genimages:=   List( mapi2[1],
			    v -> ImagesRepresentative( map1, v ) )
		    + mapi2[2];

      fi;

    fi;

    # Construct the linear mapping.
    sum:= LeftModuleHomomorphismByImagesNC(
              Source( map1 ), Range( map1 ), gens, genimages );

    # Return the sum.
    return sum;
    end );


#############################################################################
##
##  2. methods for linear mappings given by matrices
##


#############################################################################
##
#R  IsLinearMappingByMatrixDefaultRep
##
##  is another default representation of $F$-linear mappings between
##  two free left modules $V$ and $W$ where $F$ is equal to the left acting
##  domain of $V$ and of $W$.
##
##  Defining components are
##
##  `basissource' \: \\
##      basis of $V$,
##
##  `basisrange' \: \\
##      basis of $W$,
##
##  `matrix' \: \\
##      matrix over $F$, of dimensions $\dim(V)$ times $\dim(W)$.
##
##  The mapping is defined as follows.
##  The image of a vector in $V$ has coefficients
##  `Coefficients( <map>!.basissource <v> ) * <map>!.matrix'
##  w.r.t. `<map>!.basisrange'.
##
##  If one wants to compute preimages, one needs the components
##  `basisimage' \: \\
##      basis of the image of <map>,
##
##  `preimagesbasisimage' \: \\
##      preimages of the basis vectors of `basisimage',
##
##  `relations' \: \\
##      linearly independent generators for the relation space,
##      i.e., of the left null space of `<map>!.matrix'.
##      (The corresponding linear combinations of `basissource'
##      generate the kernel.)
##
##  If these components are not yet bound, they are computed by
##  `MakePreImagesInfoLinearMappingByMatrix'.
##
##  Computed images and preimages of free left modules under linear mappings
##  are always free left modules.
##  If one needs more structure (e.g., that of an algebra) for an image or
##  preimage then the linear mapping must have a special representation.
##
##  Note that the inverse general mapping of a linear mapping defined by
##  a matrix is best handled if it uses the default method,
##  since such an inverse general mapping delegates the tasks of computing
##  (pre)images to the original general mapping.
##  So the (pre)images info is computed only once.
#T  but what about sums of such mappings?
#T  better try to share info also in this case?
#T  (share a list that is filled with the info later?)
##
DeclareRepresentation(
    "IsLinearMappingByMatrixDefaultRep",
    IsAttributeStoringRep,
    [ "basissource", "basisrange", "matrix",
      "basisimage", "preimagesbasisimage", "relations" ] );

InstallTrueMethod( IsAdditiveElementWithInverse,
    IsGeneralMapping and IsLinearMappingByMatrixDefaultRep );
InstallTrueMethod( IsLeftModuleGeneralMapping,
    IsGeneralMapping and IsLinearMappingByMatrixDefaultRep );


#############################################################################
##
#M  LeftModuleHomomorphismByMatrix( <BS>, <matrix>, <BR> )
##
##  is the total and single-valued linear general mapping with <BS> a basis
##  of the source and <BR> a basis of the range, and the rows of the matrix
##  <matrix> being the coefficients vectors of the images of <BS> w.r.t.
##  <BR>.
##
InstallMethod( LeftModuleHomomorphismByMatrix,
    "for two bases of free left modules and a matrix",
    [ IsBasis, IsMatrix, IsBasis ],
    function( BS, matrix, BR )
    local S, R, map;

    S:= UnderlyingLeftModule( BS );
    R:= UnderlyingLeftModule( BR );

    # Check the arguments.
    if   Length( BS ) <> Length( matrix )  then
      Error( "<BS> and <matrix> must have the same length" );
    elif Length( BR ) <> Length( matrix[1] )  then
      Error( "<BR> and <matrix>[1] must have the same length" );
    elif LeftActingDomain( S ) <> LeftActingDomain( R ) then
      Error( "<S> and <R> must have same left acting domain" );
    fi;
#T check entries of the matrix?

    # Make the mapping.
    map:= Objectify( TypeOfDefaultGeneralMapping( S, R,
                             IsSPGeneralMapping
                         and IsSingleValued
                         and IsTotal
                         and IsLeftModuleGeneralMapping
                         and IsLinearMappingByMatrixDefaultRep ),
                     rec(
                          basissource := BS,
                          matrix      := Immutable( matrix ),
                          basisrange  := BR
                         ) );

    # return the mapping
    return map;
    end );


#############################################################################
##
#F  MakePreImagesInfoLinearMappingByMatrix( <map> )
##
##  Provide the information for computing preimages, that is, set up
##  the components `basisimage', `preimagesbasisimage', `relations'.
##
BindGlobal( "MakePreImagesInfoLinearMappingByMatrix", function( map )
    local ech,
          B;

    ech:= SemiEchelonMatTransformation( map!.matrix );
    B:= Basis( Range( map ) );
    map!.basisimage          := BasisNC( ImagesSource( map ),
                                    List( ech.vectors,
                                      x -> LinearCombination( B, x ) ) );
    map!.relations           := Immutable( ech.relations );

    map!.preimagesbasisimage := Immutable( List( ech.coeffs,
                                    row -> LinearCombination(
                                               map!.basissource, row ) ) );
end );


#############################################################################
##
#M  KernelOfAdditiveGeneralMapping( <map> ) . . . . .  for left module m.b.m.
##
InstallMethod( KernelOfAdditiveGeneralMapping,
    "for left module m.b.m.",
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( map )
    local generators, S;

    # Form the linear combinations of the basis vectors for the
    # relation space with the `basissource' of `map'.

    if not IsBound( map!.relations ) then
      MakePreImagesInfoLinearMappingByMatrix( map );
    fi;
    generators:= BasisVectors( map!.basissource );
    S:= Source( map );
    return LeftModuleByGenerators( LeftActingDomain( S ),
               List( map!.relations,
                     r -> LinearCombination( generators, r ) ),
               Zero( S ) );
    end );


#############################################################################
##
#M  IsInjective( <map> )  . . . . . . . . . . . . . .  for left module m.b.m.
##
InstallMethod( IsInjective,
    "for left module m.b.m.",
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( map )
    local generators;

    if not IsBound( map!.relations ) then
      MakePreImagesInfoLinearMappingByMatrix( map );
    fi;
    generators:= BasisVectors( map!.basissource );
    return ForAll( map!.relations,
                   r -> IsZero( LinearCombination( generators, r ) ) );
    end );


#############################################################################
##
#M  ImagesRepresentative( <map>, <elm> )  . . . . . .  for left module m.b.m.
##
InstallMethod( ImagesRepresentative,
    "for left module m.b.m., and element",
    FamSourceEqFamElm,
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep,
      IsObject ],
    function( map, elm )
    elm:= Coefficients( map!.basissource, elm );
    if elm <> fail then
      elm:= LinearCombination( map!.basisrange, elm * map!.matrix );
    fi;
    return elm;
    end );


#############################################################################
##
#M  PreImagesRepresentative( <map>, <elm> ) . . . . .  for left module m.b.m.
##
InstallMethod( PreImagesRepresentative,
    "for left module m.b.m., and element",
    FamRangeEqFamElm,
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep,
      IsObject ],
    function( map, elm )
    if not IsBound( map!.basisimage ) then
      MakePreImagesInfoLinearMappingByMatrix( map );
    fi;
    elm:= Coefficients( map!.basisimage, elm );
    if elm = fail then
      return fail;
    fi;
    return LinearCombination( map!.preimagesbasisimage, elm );
    end );


#############################################################################
##
#M  ViewObj( <map> )  . . . . . . . . . . . . . . . .  for left module m.b.m.
##
InstallMethod( ViewObj,
    "for a left module m.b.m.",
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( map )
    Print( "<linear mapping by matrix, " );
    View( UnderlyingLeftModule( map!.basissource ) );
    Print( " -> " );
    View( UnderlyingLeftModule( map!.basisrange ) );
    Print( ">" );
    end );


#############################################################################
##
#M  PrintObj( <map> ) . . . . . . . . . . . . . . . .  for left module m.b.m.
##
InstallMethod( PrintObj,
    "for a left module m.b.m.",
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( map )
    Print( "LeftModuleHomomorphismByMatrix( ",
           map!.basissource, ", ", map!.matrix, ", ",
           map!.basisrange, " )" );
    end );


#############################################################################
##
#M  NaturalHomomorphismBySubspace( <V>, <triv> )  . . . for free left modules
##
##  Return the identity mapping.
##
InstallMethod( NaturalHomomorphismBySubspace,
    "for left module and trivial left module",
    IsIdenticalObj,
    [ IsFreeLeftModule, IsFreeLeftModule and IsTrivial ],
    SUM_FLAGS, # better than everything else
    function( V, W )
    return IdentityMapping( V );
    end );


#############################################################################
##
#F  NaturalHomomorphismBySubspaceOntoFullRowSpace( <V>, <W> )
##
InstallGlobalFunction( NaturalHomomorphismBySubspaceOntoFullRowSpace,
    function( V, W )
    local F,
          Wvectors,
          mb,
          compl,
          gen,
          B,
          img,
          canbas,
          zero,
          Bimgs,
          nathom;

    # Check that the modules are finite dimensional.
    if not IsFiniteDimensional( V ) or not IsFiniteDimensional( W ) then
      TryNextMethod();
    elif not IsSubset( V, W ) then
      Error( "<W> must be contained in <V>" );
    fi;

    # If the left acting domains are different, adjust them.
    F:= LeftActingDomain( V );
    if F <> LeftActingDomain( W ) then
      F:= Intersection2( F, LeftActingDomain( W ) );
      V:= AsLeftModule( F, V );
      W:= AsLeftModule( F, W );
    fi;

    # If `V' is equal to `W', return a zero mapping.
    if Dimension( V ) = Dimension( W ) then
      return ZeroMapping( V, FullRowModule( F, 0 ) );
    fi;

    # Compute a basis of `V' through a basis of `W'.
    Wvectors:= BasisVectors( Basis( W ) );
    if IsEmpty( Wvectors ) then
      mb:= MutableBasis( F, Wvectors, Zero( W ) );
    else
      mb:= MutableBasis( F, Wvectors );
    fi;
    compl:= [];
    for gen in BasisVectors( Basis( V ) ) do
      if not IsContainedInSpan( mb, gen ) then
        Add( compl, gen );
        CloseMutableBasis( mb, gen );
      fi;
    od;
    B:= BasisNC( V, Concatenation( Wvectors, compl ) );

    # Compute the linear mapping by images.
    img:= FullRowModule( F, Length( compl ) );
    canbas:= CanonicalBasis( img );
    zero:= Zero( img );
    Bimgs:= Concatenation( List( Wvectors, v -> zero ),
                           BasisVectors( canbas ) );
    nathom:= LeftModuleHomomorphismByMatrix( B, Bimgs, canbas );
#T take a special representation for nat. hom.s,
#T (just compute coefficients, and then choose a subset ...)
    SetIsSurjective( nathom, true );

    # Enter the preimages info.
    nathom!.basisimage:= canbas;
    nathom!.preimagesbasisimage:= Immutable( compl );
#T relations are not needed if the kernel is known ?

    SetKernelOfAdditiveGeneralMapping( nathom, W );

    # Run the implications for the factor.
    UseFactorRelation( V, W, img );

    return nathom;
    end );


#############################################################################
##
#M  NaturalHomomorphismBySubspace( <V>, <W> ) . . . for two free left modules
##
##  return a left module m.b.m.
##
InstallMethod( NaturalHomomorphismBySubspace,
    "for two finite dimensional free left modules",
    IsIdenticalObj,
    [ IsFreeLeftModule, IsFreeLeftModule ],
    NaturalHomomorphismBySubspaceOntoFullRowSpace );


#############################################################################
##
#M  \*( <c>, <map> )  . . . . . . . . . . . . .  for scalar and linear m.b.m.
##
InstallMethod( \*,
    "for scalar and linear m.b.m.",
    [ IsMultiplicativeElement,
      IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( scalar, map )
    local mult,   # the multiple of `map', result
          F;      # left acting domain

    # Check the scalar.
    # (Maybe it is in fact another mapping, and we want to compose.)
    if     not IsInt( scalar )
       and not IsElmsColls( FamilyObj( scalar ),
                   FamilyObj( LeftActingDomain( Range( map ) ) ) ) then
      TryNextMethod();
    fi;

    # Construct the linear mapping (if possible).
    mult:= LeftModuleHomomorphismByMatrix(
               map!.basissource,
               scalar * map!.matrix,
               map!.basisrange );

    # If the scalar is invertible in the left acting domain of the source
    # then surjectivity and injectivity are maintained as well as the image.
    F:= LeftActingDomain( Source( map ) );
    if scalar in F and IsUnit( F, scalar ) then

      if HasIsInjective( map ) then
        SetIsInjective( mult, IsInjective( map ) );
      fi;
      if HasIsSurjective( map ) then
        SetIsSurjective( mult, IsSurjective( map ) );
      fi;

      if IsBound( map!.basisimage ) then
        scalar:= Inverse( scalar );
        mult!.basisimage          := map!.basisimage;
        mult!.preimagesbasisimage := Immutable(
            List( map!.preimagesbasisimage, v -> scalar * v ) );
        mult!.relations           := map!.relations;
      fi;

    fi;

    return mult;
    end );


#############################################################################
##
#M  AdditiveInverseOp( <map> )  . . . . . . . . . . . . . . for linear m.b.m.
##
InstallMethod( AdditiveInverseOp,
    "for linear m.b.m.",
    [ IsGeneralMapping and IsLinearMappingByMatrixDefaultRep ],
    function( map )
    local ainv;   # the additive inverse of `map', result

    # Construct the linear general mapping (if possible).
    ainv:= LeftModuleHomomorphismByMatrix(
               map!.basissource,
               AdditiveInverse( map!.matrix ),
               map!.basisrange );

    # Maintain preimages info.
    if IsBound( map!.basisimage ) then
      ainv!.basisimage          := map!.basisimage;
      ainv!.preimagesbasisimage := Immutable(
          List( map!.preimagesbasisimage, AdditiveInverse ) );
      ainv!.relations           := map!.relations;
    fi;

    # Surjectivity and injectivity are maintained.
    if HasIsInjective( map ) then
      SetIsInjective( ainv, IsInjective( map ) );
    fi;
    if HasIsSurjective( map ) then
      SetIsSurjective( ainv, IsSurjective( map ) );
    fi;

    return ainv;
    end );


#############################################################################
##
#M  CompositionMapping2( <map2>, map1> )  .  for left mod. hom. & lin. m.b.m.
##
InstallMethod( CompositionMapping2,
    "for left module hom. and linear m.b.m.",
    FamSource1EqFamRange2,
    [ IsLeftModuleHomomorphism,
      IsLeftModuleHomomorphism and IsLinearMappingByMatrixDefaultRep ],
    function( map2, map1 )
    local comp,        # composition of <map1> and <map2>, result
          BR,          # basis of the range of `map2'
          mat2;        # matrix corresponding to `map2'

    # Compute images for the generators of `map1'.
    if     IsLinearMappingByMatrixDefaultRep( map2 )
       and map1!.basisrange = map2!.basissource then

      BR   := map2!.basisrange;
      mat2 := map2!.matrix;

    else

      BR:= Range( map2 );
      if not IsFiniteDimensional( BR ) then
        TryNextMethod();
      fi;
      BR:= Basis( BR );
      mat2:= List( BasisVectors( map1!.basisrange ),
                 v -> Coefficients( BR, ImagesRepresentative( map2, v ) ) );

    fi;

    # Construct the linear mapping.
    comp:= LeftModuleHomomorphismByMatrix( map1!.basissource,
               map1!.matrix * mat2, BR );

    # Return the composition.
    return comp;
    end );


#############################################################################
##
#M  \+( <map1>, map2> ) . . . . . . . . . . . . . . . . for two linear m.b.m.
##
##  Two general mappings that respect addition can be added pointwise
##  if their images are equal and their preimages are equal.
##  The sum does also respect addition.
##
##  If both general mappings respect zero, additive inverses, scalar
##  multiplication then the sum also does.
##
BindGlobal( "SumOfMBMAndMapping", function( map1, map2 )
    local sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    fi;

    if    IsLinearMappingByMatrixDefaultRep( map2 )
       and map1!.basissource = map2!.basissource
       and map1!.basisrange  = map2!.basisrange then

      # If the bases in both mappings are the same,
      # it suffices to add the matrices.
      sum:= LeftModuleHomomorphismByMatrix(
                map1!.basissource,
                map1!.matrix + map2!.matrix,
                map1!.basisrange );

    else

      # Compute images of the generators of `map1' under `map2'.
      sum:= LeftModuleHomomorphismByMatrix(
                map1!.basissource,
                map1!.matrix
                + List( BasisVectors( map1!.basissource ),
                        v -> Coefficients( map1!.basisrange,
                                 ImagesRepresentative( map2, v ) ) ),
                map1!.basisrange );

    fi;

    # Return the sum.
    return sum;
end );

BindGlobal( "SumOfMappingAndMBM", function( map1, map2 )
    local sum;

    # Check that the linear mappings can be added.
    if    Source( map1 ) <> Source( map2 )
       or Range( map1 ) <> Range( map2 ) then
      Error( "<map1> and <map2> must have same source and range" );
    fi;

    if    IsLinearMappingByMatrixDefaultRep( map1 )
       and map1!.basissource = map2!.basissource
       and map1!.basisrange  = map2!.basisrange then

      # If the bases in both mappings are the same,
      # it suffices to add the matrices.
      sum:= LeftModuleHomomorphismByMatrix(
                map1!.basissource,
                map1!.matrix + map2!.matrix,
                map1!.basisrange );

    else

      # Compute images of the generators of `map2' under `map1'.
      sum:= LeftModuleHomomorphismByMatrix(
                map2!.basissource,
                List( BasisVectors( map2!.basissource ),
                      v -> Coefficients( map2!.basisrange,
                               ImagesRepresentative( map1, v ) ) )
                + map2!.matrix,
                map2!.basisrange );

    fi;

    # Return the sum.
    return sum;
end );

InstallOtherMethod( \+,
    "for linear m.b.m. and mapping",
    IsIdenticalObj,
    [ IsMapping and IsLinearMappingByMatrixDefaultRep,
      IsMapping ],
    SumOfMBMAndMapping );

InstallOtherMethod( \+,
    "for mapping and linear m.b.m.",
    IsIdenticalObj,
    [ IsMapping,
      IsMapping and IsLinearMappingByMatrixDefaultRep ],
    SumOfMappingAndMBM );


#############################################################################
##
#M  \+( <map1>, <map2> )  . . . . for mapping by matrix and mapping by images
##
InstallMethod( \+,
    "for linear m.b.m. and linear m.b.i.",
    IsIdenticalObj,
    [ IsMapping and IsLinearMappingByMatrixDefaultRep,
      IsMapping and IsLinearGeneralMappingByImagesDefaultRep ],
    SumOfMBMAndMapping );

InstallMethod( \+,
    "for linear m.b.i. and linear m.b.m.",
    IsIdenticalObj,
    [ IsMapping and IsLinearGeneralMappingByImagesDefaultRep,
      IsMapping and IsLinearMappingByMatrixDefaultRep ],
    SumOfMappingAndMBM );


#############################################################################
##
##  3. methods for vector spaces of linear mappings
##


#############################################################################
##
#M  NiceFreeLeftModuleInfo( <V> ) . . . . . .  for a space of linear mappings
#M  NiceVector( <V>, <v> )  . .  for space of lin. mappings, and lin. mapping
#M  UglyVector( <V>, <mat> )  . . .  for space of linear mappings, and matrix
##
InstallHandlingByNiceBasis( "IsLinearMappingsModule", rec(
    detect := function( F, gens, V, zero )
      local S, R;
      if not IsGeneralMappingCollection( V ) then
        return false;
      fi;
      gens:= AsList( gens );
      if IsEmpty( gens ) then
        S:= Source( zero );
        R:= Range(  zero );
      else
        S:= Source( gens[1] );
        R:= Range(  gens[1] );
      fi;

      # Check that the mappings have left modules as source and range.
      if    not IsLeftModule( S )
         or not IsLeftModule( R )
         or not ForAll( gens, IsMapping ) then
        return false;
      fi;

      # Check that all generators have the same source and range,
      # and that source and range are in fact left modules.
      if    ForAny( gens, map -> Source( map ) <> S )
         or ForAny( gens, map -> Range( map ) <> R ) then
        return false;
      fi;
      return true;
      end,

    NiceFreeLeftModuleInfo := function( V )
      local F, z, S, R;
      F:= LeftActingDomain( V );
      z:= Zero( V );
      S:= Source( z );
      R:= Range( z );

      # Write `S' and `R' over `F' (necessary for the nice left module).
      if LeftActingDomain( S ) <> F then
        S:= AsLeftModule( F, S );
        R:= AsLeftModule( F, R );
      fi;

      return rec( basissource := Basis( S ),
                  basisrange  := Basis( R ) );
      end,

    NiceVector := function( V, v )
      local info, M, i, c;
      info:= NiceFreeLeftModuleInfo( V );
      if     IsLinearMappingByMatrixDefaultRep( v )
         and info.basissource = v!.basissource
         and info.basisrange = v!.basisrange then
        M:= v!.matrix;
      else
        M:= [];
        for i in BasisVectors( info.basissource ) do
          c:= Coefficients( info.basisrange, ImagesRepresentative( v, i ) );
          if c = fail then
            return fail;
          fi;
          Add( M, c );
        od;
      fi;
      return M;
      end,

    UglyVector := function( V, mat )
      local info;
      info:= NiceFreeLeftModuleInfo( V );
      return LeftModuleHomomorphismByMatrix( info.basissource,
                                             mat, info.basisrange );
      end ) );


#############################################################################
##
##  4. methods for algebras of linear mappings
##


#############################################################################
##
#M  RingByGenerators( <homs> )  . . ring generated by a list of lin. mappings
##
##  If <homs> is a list of linear mappings of finite vector spaces then
##  we construct a hom algebra over the prime field.
##
InstallOtherMethod( RingByGenerators,
    "for a list of linear mappings of finite vector spaces",
    [ IsGeneralMappingCollection ],
    function( maps )
    local S;

    maps:= AsList( maps );
    if IsEmpty( maps ) then
      Error( "need at least one element" );
    fi;
    if not ForAll( maps, IsLeftModuleHomomorphism ) then
      TryNextMethod();
    fi;
    S:= Source( maps[1] );
    if     IsVectorSpace( S )
       and IsFFECollection( LeftActingDomain( S ) ) then
      return FLMLORByGenerators( GF( Characteristic( S ) ), maps );
    elif   IsVectorSpace( S )
       and IsCyclotomicCollection( LeftActingDomain( S ) ) then
      return FLMLORByGenerators( Integers, maps );
    else
      TryNextMethod();
    fi;
    end );


#############################################################################
##
#M  DefaultRingByGenerators( <maps> )  . . ring cont. a list of lin. mappings
##
##  If <maps> is a list of mappings of vector spaces then
##  we construct an algebra over the prime field.
##  (So this may differ from the result of `RingByGenerators' if the
##  characteristic is zero.)
##
InstallOtherMethod( DefaultRingByGenerators,
    "for a list of linear mappings of vector spaces",
    [ IsGeneralMappingCollection ],
    function( maps )
    local S;
    maps:= AsList( maps );
    if IsEmpty( maps ) then
      Error( "need at least one element" );
    fi;
    if not ForAll( maps, IsLeftModuleHomomorphism ) then
      TryNextMethod();
    fi;
    S:= Source( maps[1] );
    if     IsVectorSpace( S )
       and IsFFECollection( LeftActingDomain( S ) ) then
      return FLMLORByGenerators( GF( Characteristic( S ) ), maps );
    elif   IsVectorSpace( S )
       and IsCyclotomicCollection( LeftActingDomain( S ) ) then
      return FLMLORByGenerators( Rationals, maps );
    else
      TryNextMethod();
    fi;
    end );


#############################################################################
##
#M  RingWithOneByGenerators( <homs> ) . . . . . for a list of linear mappings
##
##  If <homs> is a list of linear mappings of a finite vector space then
##  we construct a hom algebra-with-one over the prime field.
##
InstallOtherMethod( RingWithOneByGenerators,
    "for a list of linear mappings of finite vector spaces",
    [ IsGeneralMappingCollection ],
    function( maps )
    local S;

    maps:= AsList( maps );
    if IsEmpty( maps ) then
      Error( "need at least one element" );
    fi;
    if not ForAll( maps, IsLeftModuleHomomorphism ) then
      TryNextMethod();
    fi;
    S:= Source( maps[1] );
    if     IsVectorSpace( S )
       and IsFFECollection( LeftActingDomain( S ) )
       and S = Range( maps[1] ) then
      return FLMLORWithOneByGenerators( GF( Characteristic( S ) ), maps );
    elif   IsVectorSpace( S )
       and IsCyclotomicCollection( LeftActingDomain( S ) ) then
      return FLMLORWithOneByGenerators( Integers, maps );
    else
      TryNextMethod();
    fi;
    end );


#############################################################################
##
#M  IsGeneratorsOfFLMLOR( <F>, <maps> )
#M  IsGeneratorsOfFLMLORWithOne( <F>, <maps> )
##
#T  check that sources and ranges coincide:
#T  if   ForAny( maps, map -> Source( map ) <> S or Range( map ) <> S ) then

#T  add implication that a FLMLOR of mappings is associative!

#T  for ideals construction, inherit the info?
#T    SetNiceFreeLeftModuleInfo( I, NiceFreeLeftModuleInfo( A ) );


#############################################################################
##
##  5. methods for full hom spaces
##


#############################################################################
##
#M  IsFullHomModule( V )  . . . . . . . . . . .  for space of linear mappings
##
InstallMethod( IsFullHomModule,
    "for space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection ],
    V -> Dimension( V ) = Dimension( UnderlyingLeftModule( NiceFreeLeftModuleInfo( V ).basissource ) )
             * Dimension( UnderlyingLeftModule( NiceFreeLeftModuleInfo( V ).basisrange ) ) );


#############################################################################
##
#M  Dimension( <M> )  . . . . . . . . . for full hom space of linear mappings
##
InstallMethod( Dimension,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    V ->   Dimension( UnderlyingLeftModule( NiceFreeLeftModuleInfo( V ).basissource ) )
         * Dimension( UnderlyingLeftModule( NiceFreeLeftModuleInfo( V ).basisrange  ) ) );


#############################################################################
##
#M  Random( <M> ) . . . . . . . . . . . for full hom space of linear mappings
##
InstallMethodWithRandomSource( Random,
    "for a random source and full hom space of linear mappings",
    [ IsRandomSource, IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    function( rs, M )
    local BS, BR;

    BR:= NiceFreeLeftModuleInfo( M );
    BS:= BR.basissource;
    BR:= BR.basisrange;

    return LeftModuleHomomorphismByMatrix( BS,
               RandomMat( rs,
                          Dimension( UnderlyingLeftModule( BS ) ),
                          Dimension( UnderlyingLeftModule( BR ) ),
                          LeftActingDomain( M ) ),
               BR );
    end );


#############################################################################
##
#M  Representative( <M> ) . . . . . . . for full hom space of linear mappings
##
##  This method is necessary for example for computing the `Zero' value of
##  <M>.  Note that <M> does in general *not* store any generators!
##
InstallMethod( Representative,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    function( M )
    local BS, BR;

    BR:= NiceFreeLeftModuleInfo( M );
    BS:= BR.basissource;
    BR:= BR.basisrange;

    return LeftModuleHomomorphismByMatrix( BS,
               NullMat( Dimension( UnderlyingLeftModule( BS ) ),
                        Dimension( UnderlyingLeftModule( BR ) ),
                        LeftActingDomain( M ) ),
               BR );
    end );


#############################################################################
##
#M  GeneratorsOfLeftModule( <V> ) . . . for full hom space of linear mappings
##
BindGlobal( "StandardGeneratorsOfFullHomModule", function( M )
    local BS, BR, R, one, m, n, zeromat, gens, i, j, gen;

    BR:= NiceFreeLeftModuleInfo( M );
    BS:= BR.basissource;
    BR:= BR.basisrange;
    R:= LeftActingDomain( M );
    one:= One( R );
    m:= Dimension( UnderlyingLeftModule( BS ) );
    n:= Dimension( UnderlyingLeftModule( BR ) );
    zeromat:= NullMat( m, n, R );
    gens:= [];
    for i in [ 1 .. m ] do
      for j in [ 1 .. n ] do
        gen:= List( zeromat, ShallowCopy );
        gen[i][j]:= one;
        Add( gens, LeftModuleHomomorphismByMatrix( BS, gen, BR ) );
      od;
    od;

    return gens;
end );

InstallMethod( GeneratorsOfLeftModule,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    StandardGeneratorsOfFullHomModule );


#############################################################################
##
#M  NiceFreeLeftModule( <M> ) . . . . . for full hom space of linear mappings
##
##  We need a special method since we decided not to store vector space
##  generators in full hom spaces;
##  note that the default methods for `NiceFreeLeftModule' are installed with
##  requirement `HasGeneratorsOfLeftModule'.
##
InstallMethod( NiceFreeLeftModule,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    function( M )
    if HasGeneratorsOfLeftModule( M ) then
      TryNextMethod();
    fi;
    GeneratorsOfLeftModule( M );
    if not HasGeneratorsOfLeftModule( M ) then
      TryNextMethod();
    fi;
    return NiceFreeLeftModule( M );
    end );


#############################################################################
##
#M  ViewObj( <M> )  . . . . . . . . . . for full hom space of linear mappings
#M  PrintObj( <M> ) . . . . . . . . . . for full hom space of linear mappings
##
BindGlobal( "ViewFullHomModule", function( M )
    local info;

    info:= NiceFreeLeftModuleInfo( M );
    if IsIdenticalObj( info.basissource, info.basisrange ) then
      Print( "End( " );
      View( LeftActingDomain( M ) );
      Print( ", " );
      View( UnderlyingLeftModule( info.basissource  ) );
      Print( " )" );
    else
      Print( "Hom( " );
      View( LeftActingDomain( M ) );
      Print( ", " );
      View( UnderlyingLeftModule( info.basissource ) );
      Print( ", " );
      View( UnderlyingLeftModule( info.basisrange  ) );
      Print( " )" );
    fi;
end );

InstallMethod( ViewObj,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ], SUM_FLAGS,
    ViewFullHomModule );

InstallMethod( PrintObj,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ], SUM_FLAGS,
    ViewFullHomModule );


#############################################################################
##
#M  \in( <v>, <V> ) . . . . . . . . . . for full hom space of linear mappings
##
InstallMethod( \in,
    "for full hom space of linear mappings",
    IsElmsColls,
    [ IsGeneralMapping,
      IsFreeLeftModule and IsGeneralMappingCollection
      and IsFullHomModule ],
    function( map, M )
    local info;
    info:= NiceFreeLeftModuleInfo( M );
    return     Source( map ) = UnderlyingLeftModule( info.basissource )
           and Range(  map ) = UnderlyingLeftModule( info.basisrange  )
           and IsLeftModuleHomomorphism( map );
    end );


#############################################################################
##
#M  IsPseudoCanonicalBasisFullHomModule( <B> )  . . . . for a full hom module
##
InstallMethod( IsPseudoCanonicalBasisFullHomModule,
    "for a basis of a full hom module",
    [ IsBasis ],
    function( B )
    local V;

    V:= UnderlyingLeftModule( B );
    if IsGeneralMappingCollection( V ) then
      if not IsFullHomModule( V ) then
        return false;
      fi;
      return BasisVectors( B ) = StandardGeneratorsOfFullHomModule( V );
    else
      Error( "do not know what a pseudo canon. basis for <V> is" );
    fi;
    end );


#############################################################################
##
#M  BasisVectors( <B> ) . . . for pseudo canonical basis of a full hom module
##
InstallMethod( BasisVectors,
    "for pseudo canonical basis of a full hom module",
    [ IsBasis and IsPseudoCanonicalBasisFullHomModule ],
    B -> StandardGeneratorsOfFullHomModule( UnderlyingLeftModule( B ) ) );


#############################################################################
##
#M  Coefficients( <B>, <m> )   for a pseudo canon. basis of a full hom module
##
InstallOtherMethod( Coefficients,
    "for pseudo canon. basis of a full hom module, and lin. mapping",
    IsCollsElms,
    [ IsBasis and IsPseudoCanonicalBasisFullHomModule, IsGeneralMapping ],
    function( B, map )
    local V, R, info;
    V:= UnderlyingLeftModule( B );
    if not IsGeneralMappingCollection( V ) then
      TryNextMethod();
    fi;
    R:= LeftActingDomain( V );
    info:= NiceFreeLeftModuleInfo( V );
    if     Source( map ) = UnderlyingLeftModule( info.basissource )
       and Range(  map ) = UnderlyingLeftModule( info.basisrange  )
       and IsLeftModuleHomomorphism( map ) then
      return Concatenation(
                  List( BasisVectors( info.basissource ),
                        v -> Coefficients( info.basisrange,
                                 ImagesRepresentative( map, v ) ) ) );
    else
      return fail;
    fi;
    end );

InstallMethod( Coefficients,
    "for pseudo can. basis of full hom module, and lin. m.b.m.",
    IsCollsElms,
    [ IsBasis and IsPseudoCanonicalBasisFullHomModule,
      IsMapping and IsLinearMappingByMatrixDefaultRep ],
    function( B, map )
    local V, R, info;
    V:= UnderlyingLeftModule( B );
    if not IsGeneralMappingCollection( V ) then
      TryNextMethod();
    fi;
    R:= LeftActingDomain( V );
    info:= NiceFreeLeftModuleInfo( V );
    if     map!.basissource = info.basissource
       and map!.basisrange  = info.basisrange then
      return Concatenation( map!.matrix );
    elif   Source( map ) = UnderlyingLeftModule( info.basissource )
       and Range(  map ) = UnderlyingLeftModule( info.basisrange  ) then
      return Concatenation(
                  List( BasisVectors( info.basissource ),
                        v -> Coefficients( info.basisrange,
                                 ImagesRepresentative( map, v ) ) ) );
    else
      return fail;
    fi;
    end );


#############################################################################
##
#M  Basis( <M> )  . . . . . . . . . . . . . . . . . . . . for full hom module
##
InstallMethod( Basis,
    "for full hom space of linear mappings",
    [ IsFreeLeftModule and IsFullHomModule ], 100,
    function( V )
    local B;
    B:= Objectify( NewType( FamilyObj( V ),
                                IsFiniteBasisDefault
                            and IsPseudoCanonicalBasisFullHomModule
                            and IsAttributeStoringRep ),
                   rec() );
    SetUnderlyingLeftModule( B, V );
    if IsFiniteDimensional( V ) then
      SetIsFinite( B, true );
      #T Check whether its is a) correct and even if so
      #T b) sometimes to expensive
      SetBasisVectors( B, GeneratorsOfLeftModule(V) );
    fi;
    return B;
    end );


#############################################################################
##
#M  Hom( <F>, <V>, <W> )
##
InstallMethod( Hom,
    "for division ring and two free left modules",
    [ IsDivisionRing, IsFreeLeftModule, IsFreeLeftModule ],
    function( F, V, W )
    local M;   # the free module record, result

    if V = W then
      return End( F, V );
    fi;
    if LeftActingDomain( V ) <> F then
      V:= AsLeftModule( F, V );
    fi;
    if LeftActingDomain( W ) <> F then
      W:= AsLeftModule( F, W );
    fi;

    M:= Objectify( NewType( CollectionsFamily( GeneralMappingsFamily(
                                ElementsFamily( FamilyObj( V ) ),
                                ElementsFamily( FamilyObj( W ) ) ) ),
                                IsFreeLeftModule
                            and IsFullHomModule
                            and IsLinearMappingsModule
                            and IsGeneralMappingCollection ),
                   rec() );

    SetLeftActingDomain( M, F );
    SetNiceFreeLeftModuleInfo( M, rec(
                                       basissource := Basis( V ),
                                       basisrange  := Basis( W ) ) );

    return M;
    end );


#############################################################################
##
#M  End( <F>, <V> ) . . . . . . . . .  for division ring and free left module
##
##  We use the generators that are also taken for full matrix FLMLORs.
##
InstallMethod( End,
    "for division ring and space of linear mappings",
    [ IsDivisionRing, IsFreeLeftModule ],
    function( F, V )
    local n,      # dimension of `V'
          i,      # loop over the rows
          gens,   # list of generators
          one,    # the identity of the field
          B,      # basis of `V'
          A;      # algebra, result

    if LeftActingDomain( V ) <> F then
      V:= AsLeftModule( F, V );
    fi;

    n:= Dimension( V );
    gens:= NullMat( n, n, F );
    gens:= [ gens, List( gens, ShallowCopy ) ];
    one:= One( F );

    # Construct the generators.
    gens[1][1][1]:= one;
    gens[2][1][n]:= one;
    for i in [ 2 .. n ] do
      gens[2][i][i-1]:= one;
    od;
    B:= Basis( V );
    gens:= List( gens, mat -> LeftModuleHomomorphismByMatrix( B, mat, B ) );

    # Construct the FLMLOR.
    A:= AlgebraWithOneByGenerators( F, gens );
    SetIsFullHomModule( A, true );
    SetNiceFreeLeftModuleInfo( A, rec(
                                       basissource := B,
                                       basisrange  := B ) );

    # Return the FLMLOR.
    return A;
    end );


#T InstallMethod( Field, true, [ IsGeneralMappingCollection ], 0, ... );
#T InstallMethod( DefaultField, true, [ IsGeneralMappingCollection ], 0, ... );
